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SPATIAL DENSITY APPROACH FOR MODELLING OF 
THE SPACE DEBRIS POPULATION 


Camilla Colombo,* Francesca Letizia* and Hugh G. Lewis* 


This article proposes a continuum density approach for space debris modelling. 
The debris population in Low Earth Orbit (LEO) is represented through its den¬ 
sity in semi-major axis, eccentricity and inclination. The time evolution of the 
density in orbital elements is modelled through the continuity equation. The per¬ 
turbing effect of aerodynamic drag is included in the divergence term, while the 
effect of fragmentation can be seen as source tern in the equation. The spatial 
density is then calculated from the orbital element density at each time. The pro¬ 
posed continuum method is used to analyse the evolution of the debris popula¬ 
tion in LEO; as initial condition the debris 2013 population is used. Then, the ef¬ 
fect of a breakup event is superimposed onto the global population of space de¬ 
bris and its effect analysed; the fragment distribution caused by the breakup up 
of satellite DMSP-F13 is considered as test case scenario. 


INTRODUCTION 

The space surrounding our planet is densely populated by an increasing number of man-made 
space debris, most of which have been generated from the break-up of operational satellites, aban¬ 
doned spacecraft or upper stages of launchers. Space debris is internationally recognised as a hazard 
to current and future space activities and space agencies are currently cooperating to identify appro¬ 
priate and sustainable space debris mitigation measures. 

The debris evolution in Low Earth Orbit (LEO) is dominated by the effects of the Earth’s oblate¬ 
ness and the atmospheric drag, which is the only natural way debris objects are removed from their 
orbits, to re-enter and bum in the atmosphere. Long-term studies of the debris environment per¬ 
form simulation of the space debris populations over 100 to 200 years to observe the effects of 
the growing space activities (e.g. launches), the uncertainty of the physical environment (e.g. at¬ 
mosphere model, changes in the Earth’s atmosphere due to the solar activities) and the spacecraft 
parameters (such as attitude, solar and drag coefficient, material deterioration), the consequences 
of fragmentation and explosion of inoperative objects and active satellites^. From the other side, 
these long-term studies aim at evaluating the efficacy of mitigation rules, such as passive disposal, 
collision avoidance manoeuvres, end-of-life guidelines, active debris removal, to reduce the risk to 
operating satellites and ensure the long term sustainability of space. 

Surveys of the existing evolution models are availablet 2-3 ^. Most of these evolutionary debris 
models f 4 ’ 5 - 6 - 7 ! use semi-analytical methods to propagate the dynamics under orbit perturbations and 
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make some assumptions on natural phenomena, the future evolution space activities, compliance 
with mitigation guidelines, and debris interaction (e.g. criteria for collision, number of fragmen¬ 
tation events per year). To ensure the analyses are robust to these uncertainties and to overcome 
the absence of a complete set of experimental data through observations several Monte Carlo runs 
are used to consider an large number of evolution scenarios l 8 - 9 - 10 - 4 !. This dramatically increases the 
computational time and limits the variety of the possible analyses. To overcome this limitation, sim¬ 
plified were proposed based on a grid discretization of the debris population in altitude bins and a 
variational approach to allow for a quick evaluation of the debris evolution fl U2 L In some cases the 
computation of the collision risk for a target spacecraft is done starting from the number of objects 
in each bin through a Poisson distribution. Some other models are instead based on a fitting process 
of the deterministic high-fidelity models t 13 ’ 14 ^. At the other end, a fully analytical model for LEO 
was proposed by Mclnnesf 15 !, borrowing the use of the continuity equation from fluid dynamics 
and planetary science. In his work the evolution of debris is described through their spatial density. 
Letizia et al. extended the continuity equation method to more than a single variable and ap¬ 
plied it to the fragment clouds generated by a single collision or fragmentation event in space. The 
knowledge of the spatial density and the distribution of relative velocities (between the cloud and a 
target spacecraft) within the cloud was used to compute the collision probability, via the kinetic gas 
theory t 17 '. In this ways, maps of collision risk can be produced in a very short computational time; 
these maps can be used for evaluating the risk on operative spacecraft. 

In this paper, we extend our previous research on the modelling of clouds through a continuum 
approach, which demonstrated to be an efficient way to propagate the density of particles in the 
space of orbital elements. We use a semi-analytical continuum density approach for debris mod¬ 
elling; the debris population in LEO is represented through its spatial density in orbital elements of 
semi-major axis, eccentricity and inclination. The time evolution of the density in orbital elements 
is modelled through the continuity equation that describes the debris flow evolution through a local 
representation via the Jacobian of the dynamics equations. With respect to existing particle-in-a-box 
approaches, where some representative objects are propagated to then rebuild the spatial density a 
posteriori, here an additional equation is added to the system dn /dt that describes the time history 
of the density of space debris in the phase space, similarly to was was done by Nazarenko f 18 ' and 
Smirnov et al. t 19 '. The proposed continuum method is validated though comparison with the actual 
debris evolution fully propagated element-wise by a semi-analytical propagator. As initial condition 
the debris population in January 2013 is used. Then, a source term is added to the continuum equa¬ 
tion, which represent a fragmentation. New fragments are thus added onto the population and their 
effect is superimposed onto the whole debris population; the case of the breakup of DMSP-F13 is 
considered. This paper will briefly describes, in the first Section, the approach developed for the 
propagation of debris fragments. The second Section will detail the application of the density based 
method to the description of the debris density evolution. The method will be applied in the third 
Section to study the evolution of the debris population in 2013 provided by the European Space 
Agency. 

DENSITY-BASED PROPAGATION FOR A CLOUD OF DEBRIS FRAGMENTS 

The propagation method ClELO (debris Cloud Evolution in Low Orbits) t 20 ' was developed to 
the aim of describing the evolution of space debris fragments resulting from breakup and collision 
in space and to assess the risk that they pose to operative spacecraft. Indeed, even in case of low 
intensity fragmentations, thousand of objects of dimension smaller than 5 cm are created. The 
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Figure 1: Schematics of the CiELO method. 


inclusion of all these objects in long-term evolutionary studies would be prohibitive in terms of 
computational time. Within our approach the fragmentation cloud is described in terms of its spatial 
density, whose evolution in time under the effect of drag is obtained by applying the continuity 
equation, following the approach proposed by Mclnnest 15 l 

The breakup is modelled through the standard NASA breakup model P 1 ’ 22 ! that gives the distribu¬ 
tion of objects in terms of their relative velocity with respect to the nominal orbital velocity where 
the fragmentation took place (which is function of the kinetic energy contribution from the event) 
and the distribution of area-to-mass. The following evolution of the fragments in LEO is dominated 
by Earth’s oblateness and atmospheric drag. In particular the effect of the Earth’s oblateness causes 
the distribution of the anomaly of the ascending node and the anomaly of the perigee of the frag¬ 
ments orbits. This phenomenon takes place over a period of time in the order of months, until the 
objects form a band around the Earth with minimum and maximum latitude approximately equal 
to the inclination where the initial fragmentation took place. For the following phase of the evo¬ 
lution, the atmospheric drag can be considered as the main perturbation and it works as a natural 
sink mechanics which removes fragments from their original orbits. In this regime, the continuous 
method can be applied to find an analytical expression which describes the time evolution of the 
spatial density. Compared to formulation by^, where the debris density is function of the radial 
distance from the Earth (r) only, the continuum method was extended to express the cloud density 
as function of semi-major axis (a) and eccentricity (e)^ 23 \ Apart of giving an insight into the evo¬ 
lution of fragments as a whole, the proposed approach drastically reduces the computational time, 
allowing the study of many fragmentation scenarios. In Fig. 1 a schematic the CiELO method is 
shown. 

Continuum approach 

The evolution of the density is obtained through the continuity equation that describes the change 
in the density of a dispersed set starting from the knowledge of the velocities of the particles. In 
particular, if n represents the fragments density, the continuity equation can be written as 

r)n 

-5- + V • f = n+ - rT (1) 


where V • f models the forces acting on the system and accounts for slow/ continuous phenomena 
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(such as orbit perturbations) and n + — h~ represents the sources and the sinks of the system, so it 
can models/ast/discontinuous events (e.g., the injection of new fragments due to launches). Once 
the initial condition for n is known, the continuity equation is used to obtain its evolution with time, 
with very low computational effort. The method was previously applied to describe the evolution 
of interplanetary dust' 24,25 !, nano-satellites constellations ' 2( 4 and high area-to-mass spacecraft' 27 '. 
The multi-dimension extension of the continuity equation Eq. 1 was fully derived in 1 16 1, following 
the approach by Gor’kavyi' 28 ^ The idea is to work in the phase space of the orbital elements by 
simply writing the divergence in rectangular coordinates; this simplifies by far the mathematical 
formulation. In the last phase of the cloud evolution, when drag is the dominant factor, the semi¬ 
major axis and the eccentricity can be chosen as phase space variables. The vector f in Eq. 1 can be 
written as a vector field with two components, respectively, the rate of variation of the semi-major 
axis a and eccentricity e caused by drag: 


f = n(a , e; t) 


f v a (a,e;t)\ 
\v e {a, e; t)) 


( 2 ) 


The expressions of the velocities were further simplified by Letizia et al. ' 16 ', to obtain an explicit 
analytical solution: 


(v a = — \J hRhP o exp ( - t, ~ j f H )/(fifl,e(a),g) ^ 

\v e = 0 

where e(a) expresses a fixed reference value of the eccentricity for each value of the semi-major 
axis. The value of e(a) was set starting from the initial distribution no(a, e). Given the expression in 
Eq. 3, the continuity equation Eq. 1 can be solved adopting the method of characteristics obtaining 
the following expression for the density: 


n(a,e;t ) = n 0 (ai,ei) 


v a (ap 

v a (a) 


(4) 


with no is the initial density at the band formation and a^, e; are two functions obtained by inverting 
the characteristics of the system at the initial timet 16 ': 


CLi (tZ, t ) 


—IT log 


exp 


a — Rjj \ 

H ) 


e i(e, t) = e. 


+ e (Rh + e(u), H) 


\IRh I 
H 


(5) 

(6) 


With Equation 4 the value of the density in the phase space at any time is known once the initial 
condition is given. As an example, Fig. 2 shows the value of the density in the phase space at the 
band formation and after 1000 days for a fragmentation at 700 km. 

Once the phase space density is known at any time t in any point of the domain, the spatial density 
can be retrieved from the phase space density by the transformations developed by Sykes' 29 ^ and 
Kessler' 3 °]. 


DENSITY PROPAGATION FOR THE WHOLE DEBRIS POPULATION 

This fully-analytical density propagation method described in the previous Section can be applied 
between 700 and 1000 kni' 16 ^, therefore it can also be employed to describe the whole LEO region. 
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a — Re [km] a — Re [km] 

Figure 2: Visualisation of cloud density (in number of fragments) following a fragmentation at 700 
km at the band formation (TB = 92 days) and after 1000 days^. 


The initial debris population at time to is known. For each object in the population we know its 
type (i.e. mission related object, payload, debris, rocket body), the area-to-mass A/M and the orbit 
condition and orbital elements. 

In time, the Earth’s oblateness causes the debris’ orbits to rotate with a precession rate that de¬ 
pends on the object’s orbital parameters and that is, therefore, different among the objects in the 
population. We can expect that, after a certain time, the right ascension of the ascending node 
and the anomaly of the perigee will be equally distributed among objects of the same kind, due to 
differences in launching time and conditions. 

As a first attempt in applying the proposed continuous technique to the global population of space 
debris, some simplifying assumptions will be made and they are justified here. The mean anomaly 
of the objects in long-term propagation studies is usually randomised, while many Monte Carlo runs 
are used to take into account differences in initial conditions, together with the uncertainties in the 
models I 8 - 9 - 10 - 4 !, in the continuum approach this is equivalent to assume that the mean anomaly of the 
object can be considered to be uniformly distributed across each orbit, therefore it can be removed as 
a variable from the continuity equation. The argument of perigee and the longitude of the ascending 
node are also randomised. Therefore, M, ui and Q can be excluded from the dependence of f in Eq. 
3. With the hypothesis of a non-rotating atmosphere, the dependence on the inclination i can also 
be removed. Under these assumption, f can be written as a vector field with two components in a 
and e as in Eq. 3. 

Eq. 4 can be now used on each point of the initial grid of the a and e domain to compute how 
the phase-space density evolve over time. Note that, in this work, we are assuming that no further 
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launches are recorded after time to, no collision among satellites are considered and no objects are 
removed from the population due to active debris removal. Each one of these terms will be added in 
a future extension of this work as the continuity equation can be also handle this cases though the 
term h + — n~ in Eq. 1. 

As said, the only perturbation on the debris population is due to the effect of drag. The effect 
will be different depending on the object’s A/M. This is tackled by dividing the considered domain 
in A/M bins. Note that in Ref. the multi-dimensional extension of the continuity equation was 
also applied to record the evolution of object with different area-to-mass; A/M, indeed can be 
added as a further phase-space variable, with a zero variation in time (i.e. the area-to-mass does not 
change over time). However, it was demonstrated that this did not result in an improvement in the 
computational time; for this reason, the binning approach is used here for the A/M. 

Fragmentation as superimposition of the effects 

Now, let’s suppose that a fragmentation takes place at a given time t,f. The NASA breakup 
model I 22 - 21 ! can be used to obtain the fragment distribution given the mass of the projectile an its 
impact velocity v c . The method previously described can be now used to compute the evolution 
of the density of the fragment cloud for any time t > tf. Therefore, at time t^ the fragments 
resulting from the fragmentation event add up to the whole debris population. From a mathematical 
point of view, this means that the new objects are added to the a — e grid and used to compute the 
new phase-space density at time tjr. The propagation is then continued to evaluate the effect of the 
fragmentation cloud on the whole debris population. 

RESULTS 

The debris population for January 2013 is used here as initial condition; this is limited to objects 
larger than 10 cm. Only objects in LEO are considered with semi-major axis a < 2000 km. Figure 
3 shows the initial distribution in semi-major axis and inclination. In Figure 4 instead, the initial 
debris distribution is shown in semi-major axis and eccentricity, distinguishing among the object 
types (MRO = Mission related Object, PL = payload, DEB = Debris, RB = Rocket Body). It can be 
noted that the objects in LEO larger than 10 cm have much lower eccentricity values than for the case 
of single breakups where the fragments have higher area-to-mass ratio. The grid considered for the 
computation is discretized with bin sizes of 20 km for semi-major axis and 0.0002 for eccentricity. 
The objects were divided in 15 bins of area to mass ratio in the rage of 0.001-13.35 m 2 /kg such 
that each bin contains the same number of objects (the difference in number of objects in each bin 
is less than 5%). 

Debris long-term evolution 

The evolution of the debris population can be computed with the method proposed, the density in 
phase space is propagated with the continuity equation, then the spatial density is calculated. Fig¬ 
ure 5 show the debris population evolution over 25 years. The results obtained with the analytical 
method (continuous line) are compared with the results obtained by a numerical method (dashed 
line) which integrates each single objects and reconstruct the spatial density though a binning ap¬ 
proach. The continuous method is able to accurately follow the debris evolution. The aerodynamic 
drag acts differently depending on the value of the area-to-mass ratio of the objects as visible from 
Figure 6. This was already noted by Mclnnest 15 !, however here the real debris population is used 
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Semi-major axis - R E [km] 

Figure 3: Initial debris distribution in LEO in semi-major axis and inclination. 
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Figure 4: Initial debris distribution in LEO in semi-major axis and eccentricity, distinguishing 
among the object types. 
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Figure 5: Debris population evolution over 25 years. Continuous line: analytical, dashed line: 
numerical 
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Figure 6: Debris population evolution over 25 years (analytical propagation). The evolution of two 
different area-to-mass bins is shown. 
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and the propagation is performed in a and e, not only in r. The method can also be applied to longer 
term studied of 100 years as shown in Figures 7 and 8. 



Figure 7: Debris population evolution over 100 years. Continuous line: analytical, dashed line: 
numerical. 



Figure 8: Debris population evolution over 100 years continuous line with analytical propagation. 


2757 








Effect of a fragmentation 

A breakup on an 800 km Sun-synchronous orbit is now considered and the feedback of the event 
on the whole debris population assessed. For the fragments distribution, the case of the fragmen¬ 
tation of DMSP-F13 on 25 February 2015 is used a test case scenario. The object had an orbit of 
perigee and apogee altitude of respectively 842 and 856 km and a mass of 830 kg. The proposed 
method can be used for a quick estimation of the future development of space debris population also 
under the effects of a fragmentation event. Figure 9 shows the spatial density of the space debris 
at the moment the fragmentation takes place and a zoom on the altitude where the fragmentation 
happens, where a distinct jump in spatial density can be recognised. The light blue line shows the 
initial population (January 2013), the darker blue line the debris evolution without fragmentation 
(in February 2015) and the red line adds up the fragmentation (in February 2015). The effect of 
the fragmentation after 20 year since January 2013 can be seen in Figure 10, a density difference of 
6.5% at the altitude of the fragmentation at time is diluted over time and in 2038 is decreased to 
3.8% (visible in the peak at altitude of around 750 km). So, over time the effect of drag level out 
the cloud peak to the envelope of the whole debris population. 
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Figure 9: Debris population + DMSP-F13 fragmentation after 2.149 years. 
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Figure 10: Debris population + DMSP-F13 fragmentation over 20 years from to. 


CONCLUSIONS 

The modelling of the contribution of small debris fragments to the collision risk requires methods 
that do not rely on the propagation of single objects; in this case, density-based models offer an 
interesting alternative. A method based on the continuity equation, previously developed to describe 
the evolution of the density of debris clouds produced by single fragmentations, is here extended to 
the study of the global debris population. The results presented show the feasibility of this approach 
for such applications with long term propagation. The accuracy of the method is demonstrated 
over long time span of 100 years so it is suitable for environment evolution studies. The effect of 
a fragmentation on the background population can be easily modelled though the superimposition 
of the effects. Future work will be devoted to complete model the sources and sinks of the debris 
population system and to measure the collision risk for spacecraft considering also background 
population. Such an approach has a potential application to perform collision risk analysis for small 
satellites. 
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